El 3 de enero de 2022, los sensores MODIS y VIIRS detectaron puntos de calor sobre el área de Macutico, en el límite entre los parques nacionales Armando Bermúdez y José del Carmen Ramírez, cordillera Central de República Dominicana. Brigadas del Ministerio de Medio Ambiente sofocaron el incendio. Puntos de calor estuvieron visibles hasta, como muy tarde, el 12 de enero.
Anomalía térmica e imagen MODIS/Terra de 3 de enero de 2022. Fuente
Usando escenas Landsat 8 (prefuego, 22/dic/2021; posfuego, 23/ene/2022), y aplicando una versión modificada del script Google Earth Engine de UN-Spider “BURN SEVERITY MAPPING USING THE NORMALIZED BURN RATIO (NBR)”, realicé un análisis de la severidad de quemado en dicho incendio usando el diferencial del índice normalizado de quema (dNBR). Adapté el script fuente a la colección más reciente de imágenes Landsat 8 (LANDSAT/LC08/C02/T1_L2).
Demo en Google Earth Engine
La web está repleta de recursos sobre esta técnica, recomiendo esta particularmente. Como resultado del análisis aplicado a Macutico, estimé unos 9.5 km2 de superficie quemada, predominantemente de baja severidad. A continuación, muestro el código de análisis y visualización de R.
library(raster)
## Loading required package: sp
library(leaflet)
dnbr <- raster('data/dNBR_macutico.tif')
Código de EarthEngine aquí.
Para aplicar a una cuenta de EarthEngine, rellenar formulario de aplicación.
dnbr_utm <- projectRaster(dnbr, crs = CRS('EPSG:32619'))
rangos <- c(
-Inf, -500, -1,
-500, -251, 1,
-251, -101, 2,
-101, 99, NA,
99, 269, 4,
269, 439, 5,
439, 659, 6,
659, 1300, 7,
1300, Inf, -1)
rclmat <- matrix(rangos, ncol = 3, byrow = TRUE)
# Reclasificar
dnbr_rcl <- reclassify(dnbr_utm, rcl = rclmat, right=NA)
# Tabla de atributos
dnbr_rcl <- ratify(dnbr_rcl)
rat <- levels(dnbr_rcl)[[1]]
# Resolucion (en km2)
resol <- res(dnbr_rcl)[1]*res(dnbr_rcl)[2]/1000000
# Etiquetas
etiq <- c('Regenerado, alta severidad', 'Regenerado, baja severidad',
'No quemado', 'Quemado, severidad baja', 'Quemado, severidad moderada-baja',
'Quemado, severidad moderada-alta', 'Quemado, severidad alta')
etiq_validas <- c(4, 5)
# Leyenda y colores
rat$leyenda <- paste(etiq[etiq_validas], round(table(dnbr_rcl[])*resol, 3), 'km cuad.')
levels(dnbr_rcl) <- rat
colores <- c("yellow2", "orange2")
leaflet() %>%
addTiles(group = 'OSM') %>%
addProviderTiles("Esri.NatGeoWorldMap", group="ESRI Mapa") %>%
addProviderTiles("Esri.WorldImagery", group="ESRI Imagen") %>%
addProviderTiles("CartoDB.Positron", group= "CartoDB") %>%
addRasterImage(dnbr_rcl, colors = colores, group = 'Severidad') %>%
addLayersControl(
overlayGroups = 'Severidad',
baseGroups = c("ESRI Imagen", "OSM", "ESRI Mapa", "CartoDB")) %>%
addLegend(
title = paste(
'Severidad. Total quemado:',
round(sum(table(dnbr_rcl[]))*resol, 2),
'km cuad.'),
pal = colorFactor(colores, rat$leyenda), values = rat$leyenda,
labels = rat$leyenda, position = 'bottomright') %>%
setView(
lat = mean(extent(dnbr)[3:4])-0.05,
lng = mean(extent(dnbr)[1:2]), zoom=11) %>%
suppressWarnings()